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Recently a local mean field theory for both eqiulibrium and transport properties of the Coulomb 
glass was proposed [A. Amir et ai, Phys. Rev. B 77, 165207 (2008); 80, 245214 (2009)]. We 
compare the predictions of this theory to the results of dynamic Monte Carlo simulations. In a 
thermal equilibrium state we compare the density of states and the occupation probabilities. We 
also study the transition rates between different states and find that the mean field rates used in the 
aforementioned papers underestimate a certain class of important transitions. We propose modified 
rates to be used in the mean field approach which take into account correlations at the minimal level 
in the sense that transitions are only to take place from an occupied to an empty site. We show 
that this modification accounts for most of the difference between the mean field and Monte Carlo 
rates. The linear response conductance is shown to exhibit the Efros-Shklovskii behavior in both the 
mean field and Monte Carlo approaches, but the local mean field method strongly underestimates 
the current at low temperatures. When using the modified rates better agreement is achieved. 

PACS numbers: 71.23.Cq, 72.15.Cz, 72.20.Ee 



I. INTRODUCTION 

Electronic states in disordered materials can be local- 
ized and at zero temperature the material can be in- 
sulating. At finite temperature, transport can still oc- 
cur trough hopping between the localized stated. The 
energy mismatch between the states must be supplied 
by emission or absorption of phonons, so the process is 
called phonon-assisted hopping. This has been studied 
for many years, one of the early major developments be- 
ing Mott's concept of variable range hopping^ (VRH) 
leading to the Mott law for the temperature dependence 
of the conductivity: 

acxexph(fo/r)i/('^+i)]. (1) 

Here T is temperature, Tq is some constant dependent 
on the material parameters while d is the dimensionality 
of the conduction problem. Mott gave only a heuris- 
tic derivation of this law based on the idea of a com- 
petition between tunneling distance and energy differ- 
ence. This was later rigorously derived using percolation 
arguments. '^I'S This puts the theory on a firm basis but 
it applies only in the case of non-interacting electrons. If 
Coulomb interactions are important, the theory has to 
be extended. The first step in this direction was the un- 
derstanding of the interaction-induced dip in the density 
of states (DOS).EHsl Based on a stability argument on 
states which are stable to all onc-particle jumps, Efros 
and ShklovskiP derived an upper bound on the density 
of states increasing as e''"^ where the excitation energy 
e is counted from the Fermi level. At finite temperatures 
the gap is smeared by thermal fluctuations.'^ 

While the arguments up to this point seem rigorous, 
one then proceeds with some type of mean field (MF) 
treatment .J^ Since the energy of a certain site depends on 
the occupancy of the other sites through the Coulomb in- 
teraction, the site energies will fiuctuate in time as jumps 



take place. Instead of following the fiuctuatig occupa- 
tion numbers, one replaces them with some average oc- 
cupancy of each site. In this way, one can repeat Mott's 
argument replacing the uniform density of states with the 
Coulomb gap form, leading to the Efros-Shklovskii (ES) 
law for conductivity;^ 

acxexph(To/T)i/2]. (2) 

While this law has been seen in many experiments, it 
is difficult to prove rigorously that a manifestly many- 
particle concept like the Coulomb gap density of states 
can be used in a single particle picture like that given by 
Mott. 

Recently there was renewed interest in these questions, 
and local mean field (LMF) theory was revisited^ and 
applied to time-dependent phenomena, see Refs. [13] for 
reviews. Having been applied to the VRH problenJSI the 
LMF theory has reproduced the ES temperature depen- 
dence ([2| of the conductance. The LMF approach differs 
form the conventional one,'^ which we refer to as ESMF, 
by two aspects. Firstly, the LMF equations solved nu- 
merically at finite temperature are intended to represent 
equilibrium state and therefore to allow for the finite- 
temperature smearing of the Coulomb gap. Contrary, 
the ESMF approach uses the ground state as a start- 
ing point and allows for finite temperatures only in the 
occupation numbers. Secondly, the LMF uses a differ- 
ent expression for the transition rate between the states, 
which neglects Coulomb correlations in the transferred 
energy that seems to be inconsistent. This issue will be 
discussed in detail in Sec. |llj Because of that, the LMF 
approach is still beyond full control, and there are still 
questions about its validity. It is therefore important to 
check the results, see if they can be trusted in some re- 
gions of parameters and in this way obtain limits of the 
validity of the mean field approximation. In this work 
we compare the MF analysis to dynamic Monte Carlo 
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simulations of the hopping process. This allows us to 
compare both equilibrium properties like the DOS and 
the occupation numbers, as well as dynamical aspects 
like the transition rates and the conductance. 

We will conclude that the LMF strongly underesti- 
mates the rate of a certain important class of transitions, 
and this leads to an underestimation of the current. We 
propose a modification of the transition rates entering 
the mean-filed scheme in the spirit of the considerations 
of Ref. El Sees. 10.1.2 and 10.2.1 (see also Refs. US]), 
which takes into account the interaction-induced correla- 
tions at the minimal level, and ensures that a transition 
will take place only from an occupied to an empty site. 
This gives much better estimates for the transition rates, 
and better agreement with the Monte Carlo simulations 
of the conductance. 

Note that recently the question of a finite temperature 
phase transition to a low temperature glass phase was 
discussed, and it was found that mean field theories^^ 
predict such a phase transi tion w hile this was not found 
in Monte Carlo simulations! ^^ * ^^ * 

The paper is organized as follows. In Sec. [IT] we recall 
the MF approach and explain our modified transition 
rates. In Sec. IIIII we outline the numerical Monte Carlo 
method we use to test the MF results. We comparw the 
properties (DOS and occupation probabilities) of equilib- 
rium states in Sec. |IV[ The transition rates are compared 
in Sec. IVl and the conductance in Sec. |Vll In Sec. IVIII 
we summarize the results. 



II. MEAN FIELD EQUATIONS 

In Ref. 12, a LMF approach was developed and then 
applied to the calculation of the conductivity in the vari- 
able range hopping regime. We give here a brief sum- 
mary of their method. We model the system as a set 
of sites, each of which can be either empty or occupied 
by one electron. For numerical convenience, the sites are 
arranged in a two dimensional square lattice with lattice 
constand d. Each site is given an energy (jji uniformly dis- 
tributed in the interval [—U, U]. The total single-particle 
(addition or subtraction) energy (SPE) of site i is then 



(3) 



Here distances rij between sites i and j are measured in 
units of d while energies are measured in units of e^/c? 
where e is the electronic charge. Background dielectric 
constant k is set to 1. rij = 0, 1 is the occupancy of site 
j. The compensating background charge v measured in 
unites of e and associated to each site is introduced to 
keep the system neutral. We have considered the case 
of half filling, so that the number of electrons is half the 
number of sites, and therefore v — 1/2. The occupation 
numbers n,, and therefore also the energies e^, fluctuate 
in time as the electrons jump between the sites. In the 



MF approximation, one replaces the fluctuating quanti- 
ties by their averages, associating to each site an average 
occupancy fi — (rii) and average energy 



T 



(4) 



The average occupancy is postulated to be given by the 
Fermi distribution at the average energy, 



f^ 



(5) 



Equations Q and ^ form a closed set, which we call 
the mean field equations. It has been showrfSl that the 
solutions of these equations give a density of states with a 
linear (in 2 dimensions) gap at the Fermi level as expected 
from the analysis of the Coulomb gap by the stability 
condition in the ground statcliil 

To calculate the conductance we must consider the 
transition rates between the different sites. If an elec- 
tron hops from site i to site j the change in the systems 
energy is 



(6) 



Here the final term arises because in the definition of 
ej it was assumed that site i was initially occupied and 
site j was initially empty. Therefore the hopping event 
creates an electron on site j and a hole on site i, which at- 
tract each other with the energy r~-^. The energy change 
Acij must be supplied by the emission or absorption of 
a phonon. The rate of such a process is given hy^ 

^(|Ae.,|/i?o)e'^''-/liV(Ae,,)|n,(l - n,) (7) 



where Eq is some energy scale characteristic of the 
electron-phonon interaction (we set it to 1 in the follow- 
ing). To is a microscopic timescale which we take as our 
unit of time, a is the localization length, which we set 
equal to the lattice constant d, and N{E) = l/(e^/-^ — 1) 
is the Planck function. The difference Ae^j > cor- 
responds to phonon absorption, while Ae^ < corre- 
sponds to phonon emission. In this case, |A^(Aeij)| ~ 
N{\Aeij\) + 1 allowing for spontaneous emission. 

In the mean field approximation the product 
\N{Aeij)\ni{l — rij) is replaced by its ensemble average, 
(|iV(Aeij)|rii(l — TT-j)), which is then decopled into a prod- 
uct of averages, \N{{Aeij))\fi[l — fj). As a result, the 
rates Tij are replaced by 



TO Eo 



\N{{Ae,,)M{l-f,). (8) 



In Ref. 141 it was argued that the energy change should 
be taken as the difference in the average energies without 
including the last term in Eq. ([6]), 



{Ae,j) = E, - E, . 



(9) 



The formal reason for omitting the self-interaction term 
is that it is necessary in order to have detailed balance in 
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equilibrium, 7^ 



7°j where the superscript indicates 
that these are the rates in equilibrium in the absence 
of an applied electric field. It is also argued that this 
is natural in a mean field approach since charge can be 
thought of as transferred continuously and this term is 
proportional to the transferred charge squared which is 
then infinitesimal. 

The above considerations seem worrying for two 
(closely related) reasons. First, the energies of the 
phonons must be calculated including this term, so it 
appears that the mean field approach will give an incor- 
rect energy balance between the electrons and the phonon 
bath. Second, there will be a number of transitions which 
are given a mean field rate very different (and much 
smaller) than the real rate. To see this point clearly, 
consider the situation where Ei < and Ej > while 
both 1^;;! > T and \Ej \ > T. Then /, w 1 and w 0. 
This means that it will often be the case that site i is 
occupied and site j empty. This is a necessary condition 
for the transition from site i to site j take place, and the 
fact that it is likely to happen is reflected in the fact that 
in (8| the factor /^(l — /j) « 1. But Ej — Ei is positive 
and larger than temperature so that the rate is strongly 
suppressed by the factor TV (l(Aejj) I) « e"'^^-^')/'^. The 
real energy change if the transition is from a configuration 
where site i is occupied and site j empty and where the 
single particle energies and ej are close to the averages 
Ei and Ej is given by Eq. Q. If the sites are so close 



that e. 



the process will be an emission process 



rather than an absorption process, and N(\Aeij\) f« 1 so 
the transition rate is exponentially larger than what is 
found using the approximation (|9|. 

It appears that for the considered configuration the ap- 
proximation (|9| is not a good one and misses an impor- 
tant physical property of the process. To improve this 
while keeping the microscopic balance we propose the 
following. We keep the mean field equations ^ and ^ 
for calculating the average occupation numbers fi. But 
when calculating the rate jij we consider the joint occu- 
pation probabilities for both sites. This follows closely 
the reasoning of Ref. [TT] (Sec 10.1.2) except that there 
all other sites were considered to be in the ground state 
configuration, whereas here we assume them to have the 
mean field average occupation. 

Let F{ni, rij) denote the probability to find the system 
with occupation numbers n.^and Uj, all other sites having 
the mean field occupation /fc. F{ni,nj) depends on the 
energy changes when adding particles to the sites i and 
j. We define 



E 



fk 



(10) 



are empty, the energy changes to the other states are 

E{ni,nj) = Ei^j^Ui + Ej(^^n.j -f n^rij/rij. (11) 
We then get 

F(n„n,) = Z^:^ie-^("-"^)/^ (12) 

where 



1 



is the partition function. The transition rate is then 
UAE, 

To 



l^|^e"2'-../a|^(Ai5 )|^(1,0). (13) 



These rates satisfy the detailed balance condition 7^" = 
7°j in the absence of an applied electric field and they 
properly take into account the energy change in the tran- 
sition. Equation (131 plus a LMF scheme takes into ac- 



count both the intra-site and inter-site correlations at 
the minimal level of ensuring that jumps only take place 
from occupied to empty sites and the finite temperature 
modification of the ground state. We refer to this theory 
as the modified local mean field (MLMF) theory. In the 
following we will compare this to Monte Carlo simula- 
tions. 



III. NUMERICAL ALGORITHM 

We used the kinetic Monte Carlo algorithm introduced 
in Refs. [Ei and [20] It consists in writing the rate j?]) as 



a product of a distance dependent part, 
and an energy dependent part. 



F^ 
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'WAe,,\/Eo)\N{Ae,,)\n,{l - n^) . 



F^ is independent of the configuration and can be pre- 
calculated and stored. Since we are working on a lattice, 
it depends only on the relative distance, and this is not 
too costly. The program first selects an electron at ran- 
dom and then a possible jump weighted by the rates F^. 
If the final site is empty, the jump is then accepted with 
probability F^. Because of the linear dependence of F-^ 
on \Aeij I it is unbounded for large negative Aeij . For this 
reason, the rate has to be limited by some maximal rate, 
which we have taken rather arbitrarily as 3T where T is 
the temperature and the rates then normalized so that 
the maximal rate has probability 1 for being accepted. 
We have checked that the exact value of the maximal rate 
is not important as long as we are not applying strong 
electric fields,!^ far beyond the linear response we will 
consider below. 



as the energy of site j not taking into account the in- 
teraction with site j (note that here we are omitting the 
background charge v since it will cancel from all energy 
differences) . Starting from the state where both the sites 



IV. EQUILIBRIUM STATES 

We can now compare several properties of the equilib- 
rium states as described by the LMF and Monte Carlo 
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methodsP2l In all comparisons we use the same sample 
for both methods, so that we are certain that all differ- 
ences are a result of the approximations of the mean field 
model. We use a sample of size 100 x 100 sites and an on- 
site disorder U = 1. In the Monte Carlo simulations we 
average over from 1 to 6 initial configurations (at high 
temperatures, thermal motion is sufhcient to wash out 
any difference between initial states). From each initial 
state we perform 3 • 10^ jumps to equilibrate the system 
and then we average over the next 5 • 10^ jumps. We 
solve the mean field equations numerically by an itera- 
tive procedure. The presented results are averages over 
100 different solutions of the mean field equations. 

Let us first study the density of single particle states 
{ei in the Monte Carlo simulations Ei in the mean field 
solutions). At low temperatures this should display the 
Coulomb gap, which according to the stability analy- 
sis of one-particle stable states should be linear in two 
dimensions.^ It was previously shown^^ that in one- 
particle stable states with U = I we get in two dimen- 
sions a Coulomb gap which is not fully linear but rather 
of the form c>c |e|^'^. With the accuracy that we are work- 
ing here we do not expect to see the departure from lin- 
ear. For the thermal equilibrium states the correspond- 
ing quadratic law in three dimensions was confirmed.^® 
The mean field equations was previously shown to give 
a linear Coulomb gajJ^^ at low temperatures. Here we 
compare the solutions of the LMF equations to Monte 
Carlo results at different temperatures (Fig. [I]) . 

At low temperatures we see that the two methods give 
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FIG. 1. The density of states of SPEs at four temperatures: 
Solid lines: Averaged Monte Carlo simulations. Dashed lines: 
Averaged mean field solutions. For Monte Caxlo simulations: 
Averaging were done over one (T = 0.1), three (T = 0.04 ), 
six (T = 0.02) and six (T = 0.01) Monte Carlo simulations;. 
Mean field solutions were averaged over 100 solutions for each 
temperature. 



similar results and close to the expected linear Coulomb 
gap. At higher temperatures the Coulomb gap is smeared 
by thermal fluctuations. We observe that the smearing 
is much more efhcicnt in the Monte Carlo simulations 
than in the mean field solutions. Thus, the mean field 
equations underestimate the smearing of the Coulomb 
gap at finite temperatures. This can also be illustrated by 
plotting the density of states at the Fermi level, g{Q), as 
function of temperature (Fig. [2]). In the LMF g{0) is close 
to a linear function, g{0) = aT. For the Monte Carlo we 
find that the dependence is slightly superlinear, but if 
fitted by a linear function we get a slope of a = 1.54 in 
reasonable agreement with the previous results, ^^^^ which 
gave a — 1.3. For the mean field we get a — 0.25, which 
is almost twice the result of Ref. where it was found 
that a = 0.15. Note however that in that work sites were 
at random. In addition, the line does not pass through 
the origin as it should, and therefore we do not believe 
that their result is numerically accurate. It used J7 = 5, 
but we have confirmed that we get the same results in 
that case so this is not the source of discrepancy. 

A self consistent equation for the energy dependence 
of the density of state^^H predicts the scaling law (in our 
units) 

g{e,T)^TJ^-2{e/T) . 

We find that our LMF results are compatible with this 
equation. However, the results of the Monte Carlo sim- 
ulations do not collapse to this law. This is also evident 
from the superlinear temperature dependence shown in 
Fig. ^ 

We can also compare the occupation numbers fi with 



■ MC, U=1 

□ MF, U=1 
MF, U=5 
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FIG. 2. Density of states at the Fermi level, g{0)- The mean 
field results are the averages over 100 {U = 1) and 50 {U = 5) 
solutions. Graphs show a near linear dependence on tempera- 
ture for mean field equations, while Monte Carlo solutions are 
near linear at higher temperatures, while showing non-linear 
behavior at low temperatures. 
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FIG. 3. Fraction of sites with intermediate occupation. 
Empty symbols: Mean field solution averaged over 100 so- 
lutions for each temperature. Filled symbols: Monte Carlo 
{rii) averaged over 5 • 10^ jumps. 



the time averaged occupation numbers {rii) in the Monte 
Carlo simulations. Figure [3] shows the fraction of sites 
which have a certain value of \ fi — 0.5| or \ {ni) — 0.5|. It 
is clear that the mean field equations underestimate the 
number of sites with intermediate occupation numbers. 



V. TRANSITION RATES 

We can compare the transition rates calculated using 
the different models. In the Monte Carlo simulations we 
simply count how often a transition takes place from a 
site of energy = to a site of energy E' — Cj. That is, 
we refer to the energies before the transition took place. 
The change in energy is then Ae = Cj — Ci — ^/rij. 

In the mean field calculations we have to take some 
care in associating the proper energies to a transition so 
that comparison to the Monte Carlo results is meaning- 
ful. The solution of the mean field equations provides 
a set of occupation numbers {fi} and energies {Ei}. 
The energy of site i before the transition is given by 
E = E^j) — Ei ~ fj/T'ij since this is the energy as- 
suming site i to be empty. The energy of site j before 
the transition is E' = Ej^^^-j + l/vij = Ej + {1 — fi)/rij 
since we know that before the transition site i is occu- 
pied. In Fig. |4] we show the average of the sum of the 
rates from energy E to E' for the three different models 
at T = 0.04. We observe clearly what we discussed in 
Sec. [nj The modified mean field rates closely follow the 
Monte Carlo results, while the original mean field rates 
are smaller for transitions crossing the Fermi level {E < Q 
and > 0). 



Monte Carlo 
0.5| 




FIG. 4. Transition rates at T=0.04. Top: Monte Carlo. 
Bottom left: LMF. Bottom right: Modified LMF. For each 
model the logarithm of the average of the sum of all rates 
from a site of energy i5 to a site of energy E' is shown. White 
indicates that the rate was smaller than the smallest colored 
rate. 



VI. CONDUCTANCE 

We also calculated the conductance in all three models. 
In the Monte Carlo simulations this was done by apply- 
ing an electric field Ef in the x-direction. This modifies 
the energy change to Ae^j = tj — Ci — l/r^j — EfAxij. 
The transition rates must then be calculated using this 
energy change but othervise the simulation is as in the 
equilibrium case. The current is measured directly as the 
transferred charge in the direction of the field. 

In the mean field calculations we find the Miller- 
Abrahams resistance^SH j^^^ _ T/'yfj and construct the 
resistance network. In this case we do not use periodic 
boundary conditions in the direction of the field. Rather, 
we set all sites on the left edge to one potential and all 
on the right to a different potential. We then solve the 
Kirchhoff equations for the network to find the current 
in each resistor. 

The temperature dependence of the conductance is 
shown in Fig.[5]for the three models. We see that all mod- 
els follow closely the Efros-Shklovskii law ([2|, except at 
high temperatures where the conductance becomes tem- 
perature independent. At high temperatures all three 
models give close to the same conductance. At lower 
temperatures the mean field approach underestimates the 
conductance. There are two reasons for this. First, the 
mean field equations underestimates the number of states 
close to the bottom of the Coulomb gap. This means 
that the number of sites which are partially occupied, 
and active in transport is underestimated. Second, the 
original mean field transition rates are too small for an 
important class of transitions as discussed in Sec. [Vj The 
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FIG. 5. Logarithm of conductivity vs. T~^^'^ for Monte 
Carlo and mean field results, as specified in legend. Solid 
lines show linear fits at low temperatures, while markers are 
the results obtained from calculations/simulations. 



modified mean field rates corrects the second problem, 
and we see that this brings the results into much closer 
agreement with the Monte Carlo simulations. Assuming 
the Efros-Shklovskii law to hold, we can extract Tq from 
linear fits to the data. We obtain (in units of / ko) 
To = 5.9 ± 0.4 (Monte Carlo), Tq = 8.0 ± 0.1 (MLMF) 
and To = 12.0 ± 0.05 (original LMF). It turns out that 
a self-consistent type of percolation approach based on 
the ESMF (see Ref. [M] and references therein) produces 
good values for Tp . In particular, for two-dimensional 
case To = 6.5 ± 1.12212SI Note that the results for To of 
the MLMF and ESMF are significantly closer to those of 
the Monte Carlo calculation than the results of the origi- 
nal LMF. The difference between the the MFLM and the 
ESMF is that the MFLM takes into account (if only ap- 
proximately, as discussed in Sec. IV I the smearing of the 
Coulomb gap at finite temperature. This increases the 
density of states close to the Fermi level, and we expect 
that it increases the conductance. This effect should be 
more pronounced at higher temperatures, and therefore 
it is natural that To is larger when this effect is included, 
and this is indeed what we see. Why the ESMF gives a 
value closer to the Monte Carlo results is not clear, and it 
would be instructive to compare the values of the conduc- 
tance, not only Tq. In any case, the difference is not very 



large, and we conclude that the smearing of the Coulomb 
gap is not of great significance for DC conductance. 



VII. SUMMARY 

We have tested the recentl y prop osed local mean field 
theory of the Coulomb glasJ^^^Sl and compared it to 
Monte Carlo simulations. We have also proposed a mod- 
ified expression for the mean field transition rates to take 
into account correlations at the most basic level that a 
jump must take place from an empty to an occupied site. 

We have found that the LMF equations underestimate 
the number of sites in the Coulomb gap at finite tem- 
peratures. They will also underestimate the number of 
sites with intermediate occupation probability, resulting 
in occupation numbers that are close to either (empty 
site) or 1 (filled site). 

The transition rates of the original mean field equa- 
tions are strongly underestimated for an important class 
of transitions, where an electron jumps from below to 
above the Fermi level, but the distance is short enough 
so that the self interaction will compensate for the in- 
crease in SPE. 

The conductance follows closely the ES law in Monte 
Carlo simulations and both the LMF calculations. How- 
ever, the LMF gives values of the current smaller than 
what is seen in the Monte Carlo simulations. The differ- 
ence is largest when using the original LMF expression 
for the rates. The MLMF rates come much closer to the 
simulation results, indicating that a major part of the 
correlations in the system can be understood in the sim- 
ple pair approximation used. More complicated correla- 
tions, involving three or more sites certainly exist, but 
seem to be of less importance. The MLMF results are 
also close to what was found using the ESMF, indicating 
that the smearing of the Coulomb gap is not important 
in determining DC conductance. This is not surprising 
since at very low temperatures, T <C Tq, the typical en- 
ergy band contributing to conductance, ~ (TTo)^/^, is 
much larger than temperature. 

This work is part of the master project of one of the au- 
thors (EB) and more details can be found in his thesis.!^ 
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